import numpy as np
import pandas as pd
from numpy import linalg as LA
import math
# Visualisation libraries
## plotly
from plotly.offline import init_notebook_mode, iplot
import plotly.graph_objs as go
import plotly.offline as py
# Graphics in retina format
%config InlineBackend.figure_format = 'retina'
import warnings
warnings.filterwarnings('ignore')
def Figs(x, y, Yexact, y_range = None):
fig = go.Figure()
fig.add_trace(go.Scatter(x = x, y = y, mode='lines', name='Approximation'))
fig.add_trace(go.Scatter(x = x, y = Yexact, mode='lines', name='Exact'))
fig['layout']['xaxis'].update(range=[x[0], x[-1]])
if y_range != None :
fig['layout']['yaxis'].update(range=y_range)
fig.update_xaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
showgrid=True, gridwidth=1, gridcolor='Lightgray',
title_text = '$x$')
fig.update_yaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
showgrid=True, gridwidth=1, gridcolor='Lightgray',
title_text = '$y$')
fig.update_xaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
fig.update_yaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
fig.update_layout(legend_orientation='v', plot_bgcolor= 'white', height= 700, width= 750)
fig.update_traces(marker_line_color= 'Black', marker_line_width=0.5, opacity=1)
fig.show()
def forward_Euler(f, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
y[0]=y_init
for i in range(N):
y[i+1]=y[i]+h*f(x[i],y[i])
return x, y
Example Conisder the initial value problem $\begin{cases}y'+2xy=xe^{-x^2},\\ y(0) = 1 \end{cases},\quad 0 \leq x \leq 1$ with exact solution $y\left(x \right) = \left(1 + \dfrac{x^{2}}{2}\right) e^{- x^{2}}$, use the forward Euler to approximate the solution of the IVP.
# f(x,y(x)):
f=lambda x, y: x*np.exp(-x**(2))-2*x*y
# f'(x,y(x)):
f1=lambda x, y: np.exp(-x**2) - 2*y - 2*(x**2)*np.exp(-x**2)
# the eact solution y(x)
y_exact= lambda x: (1+(x**2)/2)*np.exp(-x**2)
# the initial value y(0)=1
y_init=1
# the domain
x_range=[0,1]
# the step size
h=0.01
# Forward Euler
x,y =forward_Euler(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [0.5, 1.1])
it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =forward_Euler(f, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 2.4427e-02 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 2.3036e-03 | 10.60 | 1.03 |
| 2 | 1.0e-03 | 2.2904e-04 | 10.06 | 1.00 |
| 3 | 1.0e-04 | 2.2891e-05 | 10.01 | 1.00 |
Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the backward (implicit) Euler method is defiend as follows $$y_{n+1}=y_{n}+hf(x_{n+1},y_{n+1}).$$
Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the (Forward) Euler method is defiend as follows $$y_{n+1}=y_{n}+hf(x_{n},y_{n}).$$
The improved (semi-implicit) Euler method can be also is defiend as follows \begin{align} \begin{cases} z_{n}=y_{n}+hf(x_{n},y_{n})\\ y_{n+1}=y_{n}+\dfrac{h}{2}\left(f(x_{n},y_{n})+f(x_{n+1},z_{n}) \right) \end{cases} \end{align}
def Improved_forward_Euler(f, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
y[0]=y_init
for i in range(N):
z=y[i]+h*f(x[i],y[i])
y[i+1]=y[i]+(h/2)*(f(x[i],y[i])+f(x[i+1],z))
return x, y
Example: Try the previous Example, using the semi-implicit Euler method.
f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Improved_forward_Euler(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])
it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =Improved_forward_Euler(f, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 1.0421e-01 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 1.2515e-03 | 83.27 | 1.92 |
| 2 | 1.0e-03 | 1.2745e-05 | 98.20 | 1.99 |
| 3 | 1.0e-04 | 1.2768e-07 | 99.82 | 2.00 |
Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ The explicit midpoint method is given by the formula $$y_{n+1}=y_{n}+hf\left(x_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}f(x_{n},y_{n})\right).$$
def Midpoint_method(f, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
y[0]=y_init
for i in range(N):
z=y[i]+h*f(x[i],y[i])
y[i+1]=y[i]+h*f(x[i]+h/2,y[i]+(h/2)*f(x[i],y[i]))
return x, y
f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Midpoint_method(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])
it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =Midpoint_method(f, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 1.6278e-01 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 1.8020e-03 | 90.34 | 1.96 |
| 2 | 1.0e-03 | 1.8205e-05 | 98.98 | 2.00 |
| 3 | 1.0e-04 | 1.8224e-07 | 99.90 | 2.00 |
Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the implicit midpoint method is given by $$y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},{\frac {1}{2}}(y_{n}+y_{n+1})\right).$$
Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the trapezoidal rule is given by $$y_{{n+1}}=y_{n}+{\dfrac 12}h{\Big (}f(x_{n},y_{n})+f(x_{{n+1}},y_{{n+1}}){\Big )}.$$
Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the Heun's method is given by $$ y_{n+1}=y_{n}+{\dfrac {h}{4}}\left(f(x_{n},y_{n})+f(x_{n}+\dfrac{2}{3}h,y_{n}+\dfrac{2}{3}f(x_n,y_n))\right). $$
def Heun_method(f, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
y[0]=y_init
for i in range(N):
y[i+1]=y[i]+(h/4)*(f(x[i],y[i])+3*f(x[i]+(2*h/3),y[i]+(2*h/3)*f(x[i],y[i])))
return x, y
f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Heun_method(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])
it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =Heun_method(f, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 1.4418e-01 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 1.6195e-03 | 89.02 | 1.95 |
| 2 | 1.0e-03 | 1.6386e-05 | 98.83 | 1.99 |
| 3 | 1.0e-04 | 1.6406e-07 | 99.88 | 2.00 |
It is defiend as follows
$y_{n+1}=y_{n}+{\dfrac {1}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\quad \text{for }n = 0, 1, 2, 3, \ldots,$ with \begin{align} \begin{cases} k_{1}&=h\ f(t_{n},y_{n}),\\k_{2}&=h\ f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{1}}{2}}\right),\\k_{3}&=h\ f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{2}}{2}}\right),\\k_{4}&=h\ f\left(t_{n}+h,y_{n}+k_{3}\right). \end{cases} \end{align}
def runge_kutta(f, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
y[0]=y_init
for i in range(N):
k1=f(x[i],y[i]);
k2=f(x[i]+(h/2),y[i]+(h/2)*k1);
k3=f(x[i]+(h/2),y[i]+(h/2)*k2);
k4=f(x[i+1],y[i]+h*k3);
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6
return x, y
f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =runge_kutta(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])
it=3
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =runge_kutta(f, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 1.7051e-04 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 2.0347e-08 | 8380.19 | 3.92 |
| 2 | 1.0e-03 | 2.0712e-12 | 9823.76 | 3.99 |
The four-step Adams predictor-corrector method uses the four-step Adams–Bashforth and Adams-Moulton methods together:
\begin{align} z_{n+1}&=y_{n}+\dfrac{h}{24} \left(55f(x_{n},y_{n})-59f(x_{n-1},y_{n-1})+37f(x_{n-2},y_{n-2})-9f(x_{n-3},y_{n-3})\right),\\ y_{n+1}&=y_{n}+\dfrac{h}{24} \left(9f(x_{n+1},z_{n+1})+19f(x_{n},y_{n})-5f(x_{n-1},y_{n-1})+f(x_{n-2},y_{n-2})\right). \end{align}def Adams_2step(f, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
# the first two steps in Runge-Kutta
x, y=runge_kutta(f, y_init, x_range, h)
for i in range(1,N):
z=y[i]+(h/2)*(3*f(x[i],y[i])-f(x[i-1],y[i-1]))
y[i+1]=y[i]+(h/2)*(f(x[i+1],z)+f(x[i],y[i]))
return x, y
f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Adams_2step(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])
it=3
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =Adams_2step(f, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 7.1255e-02 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 1.1251e-03 | 63.33 | 1.80 |
| 2 | 1.0e-03 | 1.1780e-05 | 95.50 | 1.98 |
def Adams_4step(f, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
# the first two steps in Runge-Kutta
x, y=runge_kutta(f, y_init, x_range, h)
for i in range(3,N):
z=y[i]+(h/24)*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])-9*f(x[i-3],y[i-3]))
y[i+1]=y[i]+(h/24)*(9*f(x[i+1],z)+19*f(x[i],y[i])-5*f(x[i-1],y[i-1])+f(x[i-2],y[i-2]))
return x, y
f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Adams_4step(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])
it=3
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =Adams_4step(f, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 1.4631e-04 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 7.0867e-08 | 2064.56 | 3.31 |
| 2 | 1.0e-03 | 7.9936e-12 | 8865.40 | 3.95 |
Taylor's Method of order $n$:
\begin{align*} y(x_{i+1}) = y(x_{i})& + h f(x_{i},y(x_{i})) + \frac{h^{2}}{2} f'(x_{i},y(x_{i})) + \ldots + \frac{h^{n}}{n!}f^{(n-1)}(x_{i},y(x_{i})) \end{align*}def Taylor_Method(f,f1, y_init, x_range, h):
# Inputs
# f: the ODE y'=f(t,y(x))
# f1: the derivative of f(t,y(x))
# y0: the initial value
# x_range: interval [a,b]
# h: step size h
N = int((x_range[-1] - x_range[0])/h)
x=np.linspace(x_range[0],x_range[-1],N+1)
y=np.zeros(N+1, dtype=float)
for i in range(0,N):
y[i+1]=y[i]+h*f(x[i],y[i])+((h**2)/2)*f1(x[i],y[i])
return x, y
f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Taylor_Method(f, f1, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])
it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)
for i in range(it):
h[i]=0.1/(10**i)
for i in range(it):
x,y =Taylor_Method(f, f1, y_init, x_range, h[i])
Yexact=y_exact(x)
Norm[i]=LA.norm(y-Yexact, np.inf)
if (i>0):
Ratio[i]=Norm[i-1]/Norm[i]
LOG[i]=math.log(Ratio[i],10)
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))
| h | Norm | Ratio | Log | |
|---|---|---|---|---|
| 0 | 1.0e-01 | 3.7985e+00 | 0.00 | 0.00 |
| 1 | 1.0e-02 | 4.4081e-01 | 8.62 | 0.94 |
| 2 | 1.0e-03 | 4.4802e-02 | 9.84 | 0.99 |
| 3 | 1.0e-04 | 4.4876e-03 | 9.98 | 1.00 |